home *** CD-ROM | disk | FTP | other *** search
/ Programmer Power Tools / Programmer Power Tools.iso / progjrn / pj_6_6.arc / FCODE.ARC / LINPACKD.FOR < prev    next >
Text File  |  1988-08-09  |  20KB  |  702 lines
  1. *     a TIME function for Ryan/McFarland Fortran and Microsoft Version 4.0
  2.  
  3. *    Author:    M. Steven Baker
  4. *    Date:    September 20, 1986
  5. *
  6.        real function second()
  7.        integer*4 hh,mm,ss,hd
  8.        call gettim(hh,mm,ss,hd)
  9.        second = float(hh)*3600 + float(mm*60+ss) + float(hd)/100
  10.        end
  11.  
  12.       double precision aa(200,200),a(201,200),b(200),x(200)
  13.       double precision time(8,6),cray,ops,total,norma,normx
  14.       double precision resid,residn,eps,epslon
  15.       integer ipvt(200)
  16.       lda = 201
  17.       ldaa = 200
  18. c
  19.       n = 100
  20.       cray = .056
  21.       write(*,1)
  22.     1 format(' Please send the results of this run to:'//
  23.      $       ' Jack J. Dongarra'/
  24.      $       ' Mathematics and Computer Science Division'/
  25.      $       ' Argonne National Laboratory'/
  26.      $       ' Argonne, Illinois 60439'//
  27.      $       ' Telephone: 312-972-7246'//
  28.      $       ' ARPAnet: DONGARRA@ANL-MCS'/)
  29.       ops = (2.0d0*n**3)/3.0d0 + 2.0d0*n**2
  30. c
  31.          call matgen(a,lda,n,b,norma)
  32.          t1 = second()
  33.          call dgefa(a,lda,n,ipvt,info)
  34.          time(1,1) = second() - t1
  35.          t1 = second()
  36.          call dgesl(a,lda,n,ipvt,b,0)
  37.          time(1,2) = second() - t1
  38.          total = time(1,1) + time(1,2)
  39. c
  40. c     compute a residual to verify results.
  41. c
  42.          do 10 i = 1,n
  43.             x(i) = b(i)
  44.    10    continue
  45.          call matgen(a,lda,n,b,norma)
  46.          do 20 i = 1,n
  47.             b(i) = -b(i)
  48.    20    continue
  49.          call dmxpy(n,b,n,lda,x,a)
  50.          resid = 0.0
  51.          normx = 0.0
  52.          do 30 i = 1,n
  53.             resid = dmax1( resid, dabs(b(i)) )
  54.             normx = dmax1( normx, dabs(x(i)) )
  55.    30    continue
  56.          eps = epslon(1.0d0)
  57.          residn = resid/( n*norma*normx*eps )
  58.          write(*,40)
  59.    40    format('     norm. resid      resid           machep',
  60.      $          '         x(1)          x(n)')
  61.          write(*,50) residn,resid,eps,x(1),x(n)
  62.    50    format(1p5e16.8)
  63. c
  64.          write(*,60) n
  65.    60    format(//'    times are reported for matrices of order ',i5)
  66.          write(*,70)
  67.    70    format(6x,'dgefa',6x,'dgesl',6x,'total',5x,'mflops',7x,'unit',
  68.      $         6x,'ratio')
  69. c
  70.          time(1,3) = total
  71.          time(1,4) = ops/(1.0d6*total)
  72.          time(1,5) = 2.0d0/time(1,4)
  73.          time(1,6) = total/cray
  74.          write(*,80) lda
  75.    80    format(' times for array with leading dimension of',i4)
  76.          write(*,110) (time(1,i),i=1,6)
  77. c
  78.          call matgen(a,lda,n,b,norma)
  79.          t1 = second()
  80.          call dgefa(a,lda,n,ipvt,info)
  81.          time(2,1) = second() - t1
  82.          t1 = second()
  83.          call dgesl(a,lda,n,ipvt,b,0)
  84.          time(2,2) = second() - t1
  85.          total = time(2,1) + time(2,2)
  86.          time(2,3) = total
  87.          time(2,4) = ops/(1.0d6*total)
  88.          time(2,5) = 2.0d0/time(2,4)
  89.          time(2,6) = total/cray
  90. c
  91.          call matgen(a,lda,n,b,norma)
  92.          t1 = second()
  93.          call dgefa(a,lda,n,ipvt,info)
  94.          time(3,1) = second() - t1
  95.          t1 = second()
  96.          call dgesl(a,lda,n,ipvt,b,0)
  97.          time(3,2) = second() - t1
  98.          total = time(3,1) + time(3,2)
  99.          time(3,3) = total
  100.          time(3,4) = ops/(1.0d6*total)
  101.          time(3,5) = 2.0d0/time(3,4)
  102.          time(3,6) = total/cray
  103. c
  104.          ntimes = 10
  105.          tm2 = 0
  106.          t1 = second()
  107.          do 90 i = 1,ntimes
  108.             tm = second()
  109.             call matgen(a,lda,n,b,norma)
  110.             tm2 = tm2 + second() - tm
  111.             call dgefa(a,lda,n,ipvt,info)
  112.    90    continue
  113.          time(4,1) = (second() - t1 - tm2)/ntimes
  114.          t1 = second()
  115.          do 100 i = 1,ntimes
  116.             call dgesl(a,lda,n,ipvt,b,0)
  117.   100    continue
  118.          time(4,2) = (second() - t1)/ntimes
  119.          total = time(4,1) + time(4,2)
  120.          time(4,3) = total
  121.          time(4,4) = ops/(1.0d6*total)
  122.          time(4,5) = 2.0d0/time(4,4)
  123.          time(4,6) = total/cray
  124. c
  125.          write(*,110) (time(2,i),i=1,6)
  126.          write(*,110) (time(3,i),i=1,6)
  127.          write(*,110) (time(4,i),i=1,6)
  128.   110    format(6(1pe11.3))
  129. c
  130.          call matgen(aa,ldaa,n,b,norma)
  131.          t1 = second()
  132.          call dgefa(aa,ldaa,n,ipvt,info)
  133.          time(5,1) = second() - t1
  134.          t1 = second()
  135.          call dgesl(aa,ldaa,n,ipvt,b,0)
  136.          time(5,2) = second() - t1
  137.          total = time(5,1) + time(5,2)
  138.          time(5,3) = total
  139.          time(5,4) = ops/(1.0d6*total)
  140.          time(5,5) = 2.0d0/time(5,4)
  141.          time(5,6) = total/cray
  142. c
  143.          call matgen(aa,ldaa,n,b,norma)
  144.          t1 = second()
  145.          call dgefa(aa,ldaa,n,ipvt,info)
  146.          time(6,1) = second() - t1
  147.          t1 = second()
  148.          call dgesl(aa,ldaa,n,ipvt,b,0)
  149.          time(6,2) = second() - t1
  150.          total = time(6,1) + time(6,2)
  151.          time(6,3) = total
  152.          time(6,4) = ops/(1.0d6*total)
  153.          time(6,5) = 2.0d0/time(6,4)
  154.          time(6,6) = total/cray
  155. c
  156.          call matgen(aa,ldaa,n,b,norma)
  157.          t1 = second()
  158.          call dgefa(aa,ldaa,n,ipvt,info)
  159.          time(7,1) = second() - t1
  160.          t1 = second()
  161.          call dgesl(aa,ldaa,n,ipvt,b,0)
  162.          time(7,2) = second() - t1
  163.          total = time(7,1) + time(7,2)
  164.          time(7,3) = total
  165.          time(7,4) = ops/(1.0d6*total)
  166.          time(7,5) = 2.0d0/time(7,4)
  167.          time(7,6) = total/cray
  168. c
  169.          ntimes = 10
  170.          tm2 = 0
  171.          t1 = second()
  172.          do 120 i = 1,ntimes
  173.             tm = second()
  174.             call matgen(aa,ldaa,n,b,norma)
  175.             tm2 = tm2 + second() - tm
  176.             call dgefa(aa,ldaa,n,ipvt,info)
  177.   120    continue
  178.          time(8,1) = (second() - t1 - tm2)/ntimes
  179.          t1 = second()
  180.          do 130 i = 1,ntimes
  181.             call dgesl(aa,ldaa,n,ipvt,b,0)
  182.   130    continue
  183.          time(8,2) = (second() - t1)/ntimes
  184.          total = time(8,1) + time(8,2)
  185.          time(8,3) = total
  186.          time(8,4) = ops/(1.0d6*total)
  187.          time(8,5) = 2.0d0/time(8,4)
  188.          time(8,6) = total/cray
  189. c
  190.          write(*,140) ldaa
  191.   140    format(/' times for array with leading dimension of',i4)
  192.          write(*,110) (time(5,i),i=1,6)
  193.          write(*,110) (time(6,i),i=1,6)
  194.          write(*,110) (time(7,i),i=1,6)
  195.          write(*,110) (time(8,i),i=1,6)
  196.       stop
  197.       end
  198.       subroutine matgen(a,lda,n,b,norma)
  199.       double precision a(lda,1),b(1),norma
  200. c
  201.       init = 1325
  202.       norma = 0.0
  203.       do 30 j = 1,n
  204.          do 20 i = 1,n
  205.             init = mod(3125*init,65536)
  206.             a(i,j) = (init - 32768.0)/16384.0
  207.             norma = dmax1(a(i,j), norma)
  208.    20    continue
  209.    30 continue
  210.       do 35 i = 1,n
  211.           b(i) = 0.0
  212.    35 continue
  213.       do 50 j = 1,n
  214.          do 40 i = 1,n
  215.             b(i) = b(i) + a(i,j)
  216.    40    continue
  217.    50 continue
  218.       return
  219.       end
  220.       subroutine dgefa(a,lda,n,ipvt,info)
  221.       integer lda,n,ipvt(1),info
  222.       double precision a(lda,1)
  223. c
  224. c     dgefa factors a double precision matrix by gaussian elimination.
  225. c
  226. c     dgefa is usually called by dgeco, but it can be called
  227. c     directly with a saving in time if  rcond  is not needed.
  228. c     (time for dgeco) = (1 + 9/n)*(time for dgefa) .
  229. c
  230. c     on entry
  231. c
  232. c        a       double precision(lda, n)
  233. c                the matrix to be factored.
  234. c
  235. c        lda     integer
  236. c                the leading dimension of the array  a .
  237. c
  238. c        n       integer
  239. c                the order of the matrix  a .
  240. c
  241. c     on return
  242. c
  243. c        a       an upper triangular matrix and the multipliers
  244. c                which were used to obtain it.
  245. c                the factorization can be written  a = l*u  where
  246. c                l  is a product of permutation and unit lower
  247. c                triangular matrices and  u  is upper triangular.
  248. c
  249. c        ipvt    integer(n)
  250. c                an integer vector of pivot indices.
  251. c
  252. c        info    integer
  253. c                = 0  normal value.
  254. c                = k  if  u(k,k) .eq. 0.0 .  this is not an error
  255. c                     condition for this subroutine, but it does
  256. c                     indicate that dgesl or dgedi will divide by zero
  257. c                     if called.  use  rcond  in dgeco for a reliable
  258. c                     indication of singularity.
  259. c
  260. c     linpack. this version dated 08/14/78 .
  261. c     cleve moler, university of new mexico, argonne national lab.
  262. c
  263. c     subroutines and functions
  264. c
  265. c     blas daxpy,dscal,idamax
  266. c
  267. c     internal variables
  268. c
  269.       double precision t
  270.       integer idamax,j,k,kp1,l,nm1
  271. c
  272. c
  273. c     gaussian elimination with partial pivoting
  274. c
  275.       info = 0
  276.       nm1 = n - 1
  277.       if (nm1 .lt. 1) go to 70
  278.       do 60 k = 1, nm1
  279.          kp1 = k + 1
  280. c
  281. c        find l = pivot index
  282. c
  283.          l = idamax(n-k+1,a(k,k),1) + k - 1
  284.          ipvt(k) = l
  285. c
  286. c        zero pivot implies this column already triangularized
  287. c
  288.          if (a(l,k) .eq. 0.0d0) go to 40
  289. c
  290. c           interchange if necessary
  291. c
  292.             if (l .eq. k) go to 10
  293.                t = a(l,k)
  294.                a(l,k) = a(k,k)
  295.                a(k,k) = t
  296.    10       continue
  297. c
  298. c           compute multipliers
  299. c
  300.             t = -1.0d0/a(k,k)
  301.             call dscal(n-k,t,a(k+1,k),1)
  302. c
  303. c           row elimination with column indexing
  304. c
  305.             do 30 j = kp1, n
  306.                t = a(l,j)
  307.                if (l .eq. k) go to 20
  308.                   a(l,j) = a(k,j)
  309.                   a(k,j) = t
  310.    20          continue
  311.                call daxpy(n-k,t,a(k+1,k),1,a(k+1,j),1)
  312.    30       continue
  313.          go to 50
  314.    40    continue
  315.             info = k
  316.    50    continue
  317.    60 continue
  318.    70 continue
  319.       ipvt(n) = n
  320.       if (a(n,n) .eq. 0.0d0) info = n
  321.       return
  322.       end
  323.       subroutine dgesl(a,lda,n,ipvt,b,job)
  324.       integer lda,n,ipvt(1),job
  325.       double precision a(lda,1),b(1)
  326. c
  327. c     dgesl solves the double precision system
  328. c     a * x = b  or  trans(a) * x = b
  329. c     using the factors computed by dgeco or dgefa.
  330. c
  331. c     on entry
  332. c
  333. c        a       double precision(lda, n)
  334. c                the output from dgeco or dgefa.
  335. c
  336. c        lda     integer
  337. c                the leading dimension of the array  a .
  338. c
  339. c        n       integer
  340. c                the order of the matrix  a .
  341. c
  342. c        ipvt    integer(n)
  343. c                the pivot vector from dgeco or dgefa.
  344. c
  345. c        b       double precision(n)
  346. c                the right hand side vector.
  347. c
  348. c        job     integer
  349. c                = 0         to solve  a*x = b ,
  350. c                = nonzero   to solve  trans(a)*x = b  where
  351. c                            trans(a)  is the transpose.
  352. c
  353. c     on return
  354. c
  355. c        b       the solution vector  x .
  356. c
  357. c     error condition
  358. c
  359. c        a division by zero will occur if the input factor contains a
  360. c        zero on the diagonal.  technically this indicates singularity
  361. c        but it is often caused by improper arguments or improper
  362. c        setting of lda .  it will not occur if the subroutines are
  363. c        called correctly and if dgeco has set rcond .gt. 0.0
  364. c        or dgefa has set info .eq. 0 .
  365. c
  366. c     to compute  inverse(a) * c  where  c  is a matrix
  367. c     with  p  columns
  368. c           call dgeco(a,lda,n,ipvt,rcond,z)
  369. c           if (rcond is too small) go to ...
  370. c           do 10 j = 1, p
  371. c              call dgesl(a,lda,n,ipvt,c(1,j),0)
  372. c        10 continue
  373. c
  374. c     linpack. this version dated 08/14/78 .
  375. c     cleve moler, university of new mexico, argonne national lab.
  376. c
  377. c     subroutines and functions
  378. c
  379. c     blas daxpy,ddot
  380. c
  381. c     internal variables
  382. c
  383.       double precision ddot,t
  384.       integer k,kb,l,nm1
  385. c
  386.       nm1 = n - 1
  387.       if (job .ne. 0) go to 50
  388. c
  389. c        job = 0 , solve  a * x = b
  390. c        first solve  l*y = b
  391. c
  392.          if (nm1 .lt. 1) go to 30
  393.          do 20 k = 1, nm1
  394.             l = ipvt(k)
  395.             t = b(l)
  396.             if (l .eq. k) go to 10
  397.                b(l) = b(k)
  398.                b(k) = t
  399.    10       continue
  400.             call daxpy(n-k,t,a(k+1,k),1,b(k+1),1)
  401.    20    continue
  402.    30    continue
  403. c
  404. c        now solve  u*x = y
  405. c
  406.          do 40 kb = 1, n
  407.             k = n + 1 - kb
  408.             b(k) = b(k)/a(k,k)
  409.             t = -b(k)
  410.             call daxpy(k-1,t,a(1,k),1,b(1),1)
  411.    40    continue
  412.       go to 100
  413.    50 continue
  414. c
  415. c        job = nonzero, solve  trans(a) * x = b
  416. c        first solve  trans(u)*y = b
  417. c
  418.          do 60 k = 1, n
  419.             t = ddot(k-1,a(1,k),1,b(1),1)
  420.             b(k) = (b(k) - t)/a(k,k)
  421.    60    continue
  422. c
  423. c        now solve trans(l)*x = y
  424. c
  425.          if (nm1 .lt. 1) go to 90
  426.          do 80 kb = 1, nm1
  427.             k = n - kb
  428.             b(k) = b(k) + ddot(n-k,a(k+1,k),1,b(k+1),1)
  429.             l = ipvt(k)
  430.             if (l .eq. k) go to 70
  431.                t = b(l)
  432.                b(l) = b(k)
  433.                b(k) = t
  434.    70       continue
  435.    80    continue
  436.    90    continue
  437.   100 continue
  438.       return
  439.       end
  440.       subroutine daxpy(n,da,dx,incx,dy,incy)
  441. c
  442. c     constant times a vector plus a vector.
  443. c     jack dongarra, linpack, 3/11/78.
  444. c
  445.       double precision dx(1),dy(1),da
  446.       integer i,incx,incy,ix,iy,m,mp1,n
  447. c
  448.       if(n.le.0)return
  449.       if (da .eq. 0.0d0) return
  450.       if(incx.eq.1.and.incy.eq.1)go to 20
  451. c
  452. c        code for unequal increments or equal increments
  453. c          not equal to 1
  454. c
  455.       ix = 1
  456.       iy = 1
  457.       if(incx.lt.0)ix = (-n+1)*incx + 1
  458.       if(incy.lt.0)iy = (-n+1)*incy + 1
  459.       do 10 i = 1,n
  460.         dy(iy) = dy(iy) + da*dx(ix)
  461.         ix = ix + incx
  462.         iy = iy + incy
  463.    10 continue
  464.       return
  465. c
  466. c        code for both increments equal to 1
  467. c
  468.    20 continue
  469.       do 30 i = 1,n
  470.         dy(i) = dy(i) + da*dx(i)
  471.    30 continue
  472.       return
  473.       end
  474.       double precision function ddot(n,dx,incx,dy,incy)
  475. c
  476. c     forms the dot product of two vectors.
  477. c     jack dongarra, linpack, 3/11/78.
  478. c
  479.       double precision dx(1),dy(1),dtemp
  480.       integer i,incx,incy,ix,iy,m,mp1,n
  481. c
  482.       ddot = 0.0d0
  483.       dtemp = 0.0d0
  484.       if(n.le.0)return
  485.       if(incx.eq.1.and.incy.eq.1)go to 20
  486. c
  487. c        code for unequal increments or equal increments
  488. c          not equal to 1
  489. c
  490.       ix = 1
  491.       iy = 1
  492.       if(incx.lt.0)ix = (-n+1)*incx + 1
  493.       if(incy.lt.0)iy = (-n+1)*incy + 1
  494.       do 10 i = 1,n
  495.         dtemp = dtemp + dx(ix)*dy(iy)
  496.         ix = ix + incx
  497.         iy = iy + incy
  498.    10 continue
  499.       ddot = dtemp
  500.       return
  501. c
  502. c        code for both increments equal to 1
  503. c
  504.    20 continue
  505.       do 30 i = 1,n
  506.         dtemp = dtemp + dx(i)*dy(i)
  507.    30 continue
  508.       ddot = dtemp
  509.       return
  510.       end
  511.       subroutine  dscal(n,da,dx,incx)
  512. c
  513. c     scales a vector by a constant.
  514. c     jack dongarra, linpack, 3/11/78.
  515. c
  516.       double precision da,dx(1)
  517.       integer i,incx,m,mp1,n,nincx
  518. c
  519.       if(n.le.0)return
  520.       if(incx.eq.1)go to 20
  521. c
  522. c        code for increment not equal to 1
  523. c
  524.       nincx = n*incx
  525.       do 10 i = 1,nincx,incx
  526.         dx(i) = da*dx(i)
  527.    10 continue
  528.       return
  529. c
  530. c        code for increment equal to 1
  531. c
  532.    20 continue
  533.       do 30 i = 1,n
  534.         dx(i) = da*dx(i)
  535.    30 continue
  536.       return
  537.       end
  538.       integer function idamax(n,dx,incx)
  539. c
  540. c     finds the index of element having max. dabsolute value.
  541. c     jack dongarra, linpack, 3/11/78.
  542. c
  543.       double precision dx(1),dmax
  544.       integer i,incx,ix,n
  545. c
  546.       idamax = 0
  547.       if( n .lt. 1 ) return
  548.       idamax = 1
  549.       if(n.eq.1)return
  550.       if(incx.eq.1)go to 20
  551. c
  552. c        code for increment not equal to 1
  553. c
  554.       ix = 1
  555.       dmax = dabs(dx(1))
  556.       ix = ix + incx
  557.       do 10 i = 2,n
  558.          if(dabs(dx(ix)).le.dmax) go to 5
  559.          idamax = i
  560.          dmax = dabs(dx(ix))
  561.     5    ix = ix + incx
  562.    10 continue
  563.       return
  564. c
  565. c        code for increment equal to 1
  566. c
  567.    20 dmax = dabs(dx(1))
  568.       do 30 i = 2,n
  569.          if(dabs(dx(i)).le.dmax) go to 30
  570.          idamax = i
  571.          dmax = dabs(dx(i))
  572.    30 continue
  573.       return
  574.       end
  575.       double precision function epslon (x)
  576.       double precision x
  577. c
  578. c     estimate unit roundoff in quantities of size x.
  579. c
  580.       double precision a,b,c,eps
  581. c
  582. c     this program should function properly on all systems
  583. c     satisfying the following two assumptions,
  584. c        1.  the base used in representing dfloating point
  585. c            numbers is not a power of three.
  586. c        2.  the quantity  a  in statement 10 is represented to 
  587. c            the accuracy used in dfloating point variables
  588. c            that are stored in memory.
  589. c     the statement number 10 and the go to 10 are intended to
  590. c     force optimizing compilers to generate code satisfying 
  591. c     assumption 2.
  592. c     under these assumptions, it should be true that,
  593. c            a  is not exactly equal to four-thirds,
  594. c            b  has a zero for its last bit or digit,
  595. c            c  is not exactly equal to one,
  596. c            eps  measures the separation of 1.0 from
  597. c                 the next larger dfloating point number.
  598. c     the developers of eispack would appreciate being informed
  599. c     about any systems where these assumptions do not hold.
  600. c
  601. c     *****************************************************************
  602. c     this routine is one of the auxiliary routines used by eispack iii
  603. c     to avoid machine dependencies.
  604. c     *****************************************************************
  605. c
  606. c     this version dated 4/6/83.
  607. c
  608.       a = 4.0d0/3.0d0
  609.    10 b = a - 1.0d0
  610.       c = b + b + b
  611.       eps = dabs(c-1.0d0)
  612.       if (eps .eq. 0.0d0) go to 10
  613.       epslon = eps*dabs(x)
  614.       return
  615.       end
  616.       subroutine dmxpy (n1, y, n2, ldm, x, m)
  617.       double precision y(*), x(*), m(ldm,*)
  618. c
  619. c   purpose:
  620. c     multiply matrix m times vector x and add the result to vector y.
  621. c
  622. c   parameters:
  623. c
  624. c     n1 integer, number of elements in vector y, and number of rows in
  625. c         matrix m
  626. c
  627. c     y double precision(n1), vector of length n1 to which is added 
  628. c         the product m*x
  629. c
  630. c     n2 integer, number of elements in vector x, and number of columns
  631. c         in matrix m
  632. c
  633. c     ldm integer, leading dimension of array m
  634. c
  635. c     x double precision(n2), vector of length n2
  636. c
  637. c     m double precision(ldm,n2), matrix of n1 rows and n2 columns
  638. c
  639. c ----------------------------------------------------------------------
  640. c
  641. c   cleanup odd vector
  642. c
  643.       j = mod(n2,2)
  644.       if (j .ge. 1) then
  645.          do 10 i = 1, n1
  646.             y(i) = (y(i)) + x(j)*m(i,j)
  647.    10    continue
  648.       endif
  649. c
  650. c   cleanup odd group of two vectors
  651. c
  652.       j = mod(n2,4)
  653.       if (j .ge. 2) then
  654.          do 20 i = 1, n1
  655.             y(i) = ( (y(i))
  656.      $             + x(j-1)*m(i,j-1)) + x(j)*m(i,j)
  657.    20    continue
  658.       endif
  659. c
  660. c   cleanup odd group of four vectors
  661. c
  662.       j = mod(n2,8)
  663.       if (j .ge. 4) then
  664.          do 30 i = 1, n1
  665.             y(i) = ((( (y(i))
  666.      $             + x(j-3)*m(i,j-3)) + x(j-2)*m(i,j-2))
  667.      $             + x(j-1)*m(i,j-1)) + x(j)  *m(i,j)
  668.    30    continue
  669.       endif
  670. c
  671. c   cleanup odd group of eight vectors
  672. c
  673.       j = mod(n2,16)
  674.       if (j .ge. 8) then
  675.          do 40 i = 1, n1
  676.             y(i) = ((((((( (y(i))
  677.      $             + x(j-7)*m(i,j-7)) + x(j-6)*m(i,j-6))
  678.      $             + x(j-5)*m(i,j-5)) + x(j-4)*m(i,j-4))
  679.      $             + x(j-3)*m(i,j-3)) + x(j-2)*m(i,j-2))
  680.      $             + x(j-1)*m(i,j-1)) + x(j)  *m(i,j)
  681.    40    continue
  682.       endif
  683. c
  684. c   main loop - groups of sixteen vectors
  685. c
  686.       jmin = j+16
  687.       do 60 j = jmin, n2, 16
  688.          do 50 i = 1, n1
  689.             y(i) = ((((((((((((((( (y(i))
  690.      $             + x(j-15)*m(i,j-15)) + x(j-14)*m(i,j-14))
  691.      $             + x(j-13)*m(i,j-13)) + x(j-12)*m(i,j-12))
  692.      $             + x(j-11)*m(i,j-11)) + x(j-10)*m(i,j-10))
  693.      $             + x(j- 9)*m(i,j- 9)) + x(j- 8)*m(i,j- 8))
  694.      $             + x(j- 7)*m(i,j- 7)) + x(j- 6)*m(i,j- 6))
  695.      $             + x(j- 5)*m(i,j- 5)) + x(j- 4)*m(i,j- 4))
  696.      $             + x(j- 3)*m(i,j- 3)) + x(j- 2)*m(i,j- 2))
  697.      $             + x(j- 1)*m(i,j- 1)) + x(j)   *m(i,j)
  698.    50    continue
  699.    60 continue
  700.       return
  701.       end
  702.